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Abstract 

During each ovarian cycle, only a definite number of follicles ovulate, while the others undergo a 
degeneration process called atresia. We have designed a multi-scale mathematical model where ovulation 
and atresia result from a hormonal controlled selection process. A 2D-conservation law describes the age 
and maturity structuration of the follicular cell population. In this paper, we focus on the operating 
mode of the control, through the study of the characteristics associated with the conservation law. We 
describe in particular the set of microscopic initial conditions leading to the macroscopic phenomenon of 
either ovulation or atresia, in the framework of backwards reachable sets theory. 

Keywords : biomathematics, conservation laws, method of characteristics, control theory, backwards 
reachable sets 

1 Introduction 

The development of ovarian follicles is a crucial process for reproduction in mammals, as its biological meaning 
is to free fcrtilizablc oocyte(s) at the time of ovulation. A better understanding of follicular development is 
both a clinical and zootechnical challenge; it is required to improve the control of anovulatory infertility in 
women, as well as ovulation rate and ovarian cycle chronology in domestic species. 

Within all the developing follicles, very few actually reach the ovulatory size; most of them undergo a 
degeneration process, known as atresia pQ. The ovulation rate (number of ovulatory follicles per cycle) 
results from an FSH (Follicle Stimulating Hormone)-dependent follicle selection process. FSH acts on the cells 
surrounding the oocyte, the granulosa cells, and controls both their commitment towards cither proliferation, 
differentiation or apoptosis and their ability to secrete hormones such as estradiol. The whole estradiol output 
from the ovaries is responsible for exerting a negative feedback on FSH release. Following the subsequent 
fall in plasmatic FSH levels, most of the follicles undergo atresia and only the ovulatory ones survive in the 
FSH-poor environment. 

We have proposed a mathematical model, using both multi-scale modeling and control theory concepts, 
to describe the follicle selection process [2]. For each follicle, the cell population dynamics is ruled by a 
conservation law with variable coefficients, which describes the changes in age and maturity of the granulosa 
cell density. A control term, representing FSH signal, intervenes both in the velocity and loss terms of the 
conservation law. Two acting controls are distinguished: a global control resulting from the ovarian feedback 
and corresponding to FSH plasmatic levels, and a local control, specific to each follicle, accounting for the 
modulation in FSH bioavaibility due to follicular vascularization. Both ovulation triggering and follicular 
ovulation depend on the reaching of a target. 

In this paper, we aim at using control theory to characterize follicular trajectories and the control laws leading 
to ovulation or atresia. The macroscopic phenomenon of follicular ovulation is considered as a reachability 
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problem for the microscopic characteristics associated with the conservation law. The problem is enounced 
in details and solved in open-loop. Follicles are assumed to ovulate if their state variables, (age, maturity and 
cell density) reach a given target set. Backwards reachable sets theory is used to define the initial conditions 
compatible with latter ovulation, for a given control panel. In section 2, the conservation laws describing 
the model for follicle selection and their characteristics are presented. In section 3, the target corresponding 
to follicular ovulation is defined. Reachable sets theory is used to solve the corresponding control problem 
and simulation results arc discussed. Section 4 is devoted to the discussion and perspectives. 



2 Follicle selection model 

2.1 Controlled conservation laws for granulosa cells 

Following [2] , cells are characterized by their positions within or outside the cell cycle and by their sensitivities 
to FSH. This leads to distinguish 3 cellular phases within the granulosa cell population. Phases 1 and 2 
correspond to the proliferation phases (describing respectively the Gl phase and S to M phases of the cell 
cycle), and phase 3 corresponds to the differentiation phase, after cells have exited the cell cycle (see Figure 
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Figure 1: Cell flow chart. The cell cycle consists of the cyclic G1-SM-G1 pathway. When it enters Gl from 
SM, a mother cell gives birth to two daughter cells. Cell differentiation corresponds to the one-way flow 
from Gl into D. Entry into apoptosis arises from the Gl and D phases. The Gl and D phases arc under the 
control of FSH signal. 

More precisely, the position of a cell at time t G [0, T] (T > 0) is defined by its age and maturity. The cell age 
a is a marker of progression within the cell cycle (in phases 1 and 2) and evolves as time t outside the cycle 
(in phase 3). The duration of phase 1 is a x > and the total cycle duration is a 2 (so that phase 2 duration 
is a 2 — ai > 0). The maturity marker 7 is used to sort the cycling and non-cycling cells by comparison to a 
threshold 7^ and to characterize the cell vulnerability towards apoptosis (programmed cell death). Phases 
1, 2, 3 correspond respectively to ranges in the values of (a, 7) in the following open sets of the age-maturity 
plane: (see Figure [2]): 

=}ka 2 , ka 2 + ai[x]0,7 s [, Q 2 ,fc =}ka 2 + 01, (k + l)a 2 [x]0, y„[, £l 3 , fe =]0, (k + l)a 2 [x}j s ,j max [ 

for integers k = 0, 1, We will also use the notation Qj,k = ^j.fc><]0, T[ for j = 1, 2, 3. 

The cell population in a follicle, /, is represented by cell density functions, -(a, 7,£), defined on each 
cellular phase Qj.k, j = 1, 2, 3, k = 0, 1, . . . as solutions of the following conservation laws [2]: 

dd> k t ■ dgf(uf)6 k f , dhfh,Uf)6 k f- 

J^L + //M ■ + nh U1M = -A( 7 ,U)<t>% mQ jtk , j = 1,2,3, * = 0,1,... (1) 
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Figure 2: Cellular phases on the age-maturity plane. The cell cycle is represented by the pathway fii,fc, 
il2,fc... and the differentiation phase by fi3,./v- 

In phases 1 and 3, both a global control, {/, and a local control, it/, act on the velocities of aging (gf 
function, locally controlled in phase 1) and maturation (hf function, locally controlled in phases 1 and 3) as 
well as on the loss term (A apoptosis rate, globally controlled in phases 1 and 3). Phase 2 is uncontrolled and 
corresponds to completion of mitosis after a pure delay in age, a 2 — a\ (no cell can leave the cycle here). We 
have the following expressions for those functions, for k = 0, 1, . . . (all parameters are real positive numbers 
defined in Table Q}: 



.9/(7,"/) = T gf (l -g l uj s _{ 1 ){\ - u f )) in Q j>k , j = 1,3 (2) 

h f{l, u f) = T hf{-1 2 + (ci7 + c 2 )(l - w s (7)cxp(-u//u))) inQ^fc, j = l,3 (3) 

A( 7 , U) - w A ( 7 )(l - U) in Q jik , j = 1, 3 (4) 

9f(u f ) = T gf , hf(-Y,u f ) = \{i,U) = mtt 2 , k (5) 



The lo s _ and tu s+ functions are smooth switches between and 1 in a strip defined for all a, < r y s _ < 7 < 7 S+ 
containing 7^ (we use the notations \x\ + = max(x,0) and = max(— x,0)): 

w s _ (7) = exp (-— — -M^) , w s+ (7) = exp (-— — ^4^) , w s (7) = (7) + u s+ (7). 
V 17-7*1 / V |7-7«I / 

The factor oj s _ in @, switches off the control of gf in phase 3: gf — r g f in ^3.^-, and switches it on for (a, 7) 
strictly in phase 1, where gf(uf) — T g f(l — .91 + fliM/)) for 7 < 7 S _ . 

In the same manner, lj s in ([3]), switches off the control of hf for 7 = j s , so that, in a neighbourhood of 7^, 
hf > for all t and uf. the boundary 7 = 7 S is inward for f23 i00 and outward for each Ox^. In particular, a 
cell very close to exit the cycle (in phase 1) is committed to go into phase 3. The exiting flux is controlled 
by Uf inside phase 1 before the cell maturity reaches 7^: the sign of hf can be changed by the control only 
below 7 S _ . 

In complement to this local control (by Uf) of the exit flux from the cell cycle, a global control is exerted 
(by U, U < 1) on the cell vulnerability towards apoptosis, in a zone surrounding 7 S , through the function 

w A (7) = i^cxp (- (^) 2 ) in ©. 

Each boundary segment of each (rectangular) O^fc, j = 1,3, k = 0,1,... is either inward or outward, 
independently on time t or controls Uf, U , so that we can consider the following boundary conditions (the 
use of the trace values and the wellposedncss of this model will be discussed later): 

Boundary conditions for (ffj x , t (=]0, T[ : 

gf^Uf^U^t) = {^ 1 i*^' 7 ' t) ' far *- 1 ' for 7 E]0, 7s [. (6) 
<^ ;1 (a,0,i) = 0, for a e]fca 2 , ka 2 + Oi[ (7) 
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Boundary conditions for 0j 3 , N > 1, t g]0, T[ : 

</>^ 3 (0,7,t) = 0, for 7 e]7 s ,7 mQ x[ (8) 

^ 8 («,7«*) = ( fJ,i(^'.;|> oe f^+ fl '[ , forfe = 0,...iV-l. (9) 
/,3V ' ' ' ' \ 0, for a e [ka 2 + a u (k + l)a 2 [ w 

Initial conditions: 

^(o,7,0) = 4>% (a, 7) | %fe , J = 1,3. (10) 

In (|6|), the delay ai — a\ and the doubling of the flux due to mitosis have been explicitly taken into account, 
so that it is not necessary anymore to consider phase 2 and (fJj 2 in the model. The cell density of interest in 

N 
f 

<$ = <t>) A on Q 1)k , fc = 0, . . . TV — 1, 4>f = <f>f j3 on Q 3 ,jv 



a follicle, ^ , is simply 



In practice, we will fix the value of N to a large enough integer, and use the notations (/>/ and a max = Na,2- 
<j)f can be computed by solving the series of problems: 

1. Pi,k- ®, @, ©, (HHD on Qi, fc , for fc = 0, . . . /V - 1. (11) 

2. P 3 ,n ■ HI), ©, ©, (HHD on Q 3 ,jv, when the P ljk arc solved . (12) 

It is worth noticing that those problems can be solved in sequence when the traces of the solutions of P\ t k 
on their outward boundaries (left and upper segments) are well defined. This is a useful property for the 
efficiency of computations and the backward reachability technique. 

The feedback exerted by the ovaries on the secretion of the hormonal control FSH defines a closed-loop 
system (cf. Figure [3]). The global control U can be interpreted as FSH plasmatic level. The local control Uf 
represents intra-follicular bioavailablc FSH levels. It is given as a proportion of U. 
Define the maturity operator M as: 

M(<p)(t) = I 1 ™* f maX 7 ^(a,7,i)darf7 (13) 
Jo Jo 

The global maturities M(<pf) on the follicular scale, and M{S~]<f)f) on the ovarian scale, will be used to 

/ 

define the two-scale feedback control: 

U = S T (M(Y^M) + U 
f 

u f = b fTf (M((t> f )).U (14) 

where S(fi) = U s + 1 ~ ?' -, S T (fi)(t) = S(fx(t - r)) 

1 + cxp (c(/x — m)) 

6/(/i) = 6l+ l + exp 1 (6 2 (t-M)) ' b t"M®= b 'W-Tf)) 

The decreasing sigmoid function S accounts for the ovarian negative feedback on FSH ; Ua(t) is a potential 
exogenous entry. The increasing sigmoid function bf makes Uf increase with the maturity of follicle /, up to 
its reaching FSH plasmatic values. The delay r is introduced to take into account the time needed for FSH 
signal to reach the ovaries; and the local delay Tf, that can be neglected in practice {jf <C t) is mainly here 
to facilitate the resolution on finite time intervals. The delayed versions of S and bf are denoted S T and 6/ T/ 
respectively. The parameters are chosen such that U and bf arc normalized: < b\, U s < 1, so that: 

< u f < U < 1 + U . (15) 
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Functions and parameters 


Definition 


Nominal value 




Global feedback control 






minimal plasmatic FSH value 


0.5 


c 


slope parameter of the sigmoid function 


0.1 




abscissa of the inflexion point of the sigmoid function 


50 


r 


FSH delay from the pituitary gland to the ovaries 


0.01 


6/[M(^/)] 


Local feedback gain 




61 


basal level 


0.054 


f>2 


exponential rate 


0.3 


63 


scaling parameter 


27 


T / 


delay to modify local vascularization 


0.01 




Aging velocity 






time scale parameter 


1 


91 


control gain in phase 1 


0.5 




Maturation velocity 




T hf 


time scale parameter 


0.07 


Cl 


slope parameter 


11.892 


C2 


origin ordinate 


2.288 




scaling parameter 


0.133 




Global feedback gain 




A" 


amplification constant 


3 


7 


scaling parameter 


0.2 


ai 


cellular age at the end of phase 1 


1 


(1-2 


cellular age at the end of phase 2 


2 


N 


maximum number of cell cycles 


8 


Is 


maturity threshold for cell cycle exit 


3 




maturity threshold for switching off control in phase 1 


2.99 


7s + 


maturity threshold for switching on control in phase 3 


3.01 


/max 


maximum maturity 


15 


M s 


ovarian threshold for ovulation triggering 


75 


M sl 


follicular threshold for ovulation ability 


40 



Table 1: Main model functions and parameters 



Ovulation is triggered when estradiol levels reach a threshold value M s . As estradiol secretion is related to 
maturity (see [2]), the ovulation time T s is defined as: 

T s = min{T | 4> f ){T) = M s } (16) 

/ 

The follicles are then sorted according to their individual maturity. The ovulatory follicles are those whose 
maturity at time T has overpassed a threshold M s \ such as M s \ < M s . The ovulation rate is computed as: 

N S , S1 = Card{/ I M{4> f )(T s ) > M sl } (17) 

For instance, Figure [4] shows, on the age-maturity domain, the cell density of either an ovulatory follicle (left 
panel) or an atretic one (right panel) at ovulation time. The cell cycle is implemented on a periodic domain, 
where the age a is reset after the cells go through mitosis. The granulosa cell repartition in the ovulatory 
follicle is characterized by a roughly older age range, and a higher maturity and cell density ranges than in 
the atretic follicle. 
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Figure 3: Ovarian system: the multi-scale model in closed-loop. The global feedback control U is modulated 
on the follicle scale (Fol i), by the follicular maturity, M((j>i), to define the local control Ufi- 



2.2 Wellposedness of the model 

We have obtained a compact expression of the model by using a density representation of the cell population 
instead of a discrete model with a huge number of cells. Even if reasoning about controlled trajectories of 
individual cells will be useful in the sequel, the global numerical simulations are based on this compact form. 
It is then useful to discuss the wellposedness of the model in order to check, to some extent, its consistency. 
Due to the delays in the feedback control laws, we can consider the it/ and U control terms as given functions 
of time on each interval of size r/. The velocities and loss term are then given functions 5/(7, t), ft/ (7, t) and 
A(7,t). In this case, we show below that there exists a unique solution to the series of problems (TTTj) . (|12p . 
In each domain Qi t k, Q3,n, the velocity field (<?/,/i/) and the loss rate A can be considered, without loss 
of expressiveness for the model, as smooth functions. Yet, initial and boundary values may be irregular, so 
that we need to check the existence of weak solutions. 

As mentioned below, it is sufficient to solve the sequence of Pi,*; and then P3,jv to define </>/. Since those 
problems have a common structure, we just need to study the solution of P\.k and its traces on the outward 
boundary used to define both Pi,fc+i and P$,n- 

We first use a change of variables to transform (TT]) into a conservation equation without a loss term. Let 
(f>j 1 = <f> H -y exp (7(7, £)) defined on Qi,u- Let dQ^ k and dQ^ k denote the inward and outward boundaries of 
the domain Qi ; k, and let l("f,t) such as 

^ + M7,t)^ = -A(7,t), 1(7,0) =0 (18) 

Here, we can assume that hf and A (extended to all 7 € K) are C 1 functions, so that, according to standard 
theory, there exists a unique generalized solution I of the linear transport equation (|18p , that is Lipschitz for 

alH > 0. Now 4> k j 1 verifies 



$1 



dt 

where vt is the (smooth) velocity vector: Vf = 



div(v f (j) L fl ) = (19) 



.9/(7,*) 
h f (j,t) 

Non homogeneous boundary conditions on dQ^ k , similar to ©, J7]) complete (fTi?)) . To solve such a problem, 
with L°° or L 1 boundary and initial values, extensions to the classical Kruzkov entropy solutions ([3]) to 



ti 



Figure 4: Cell density of an ovulatory follicle (left) and an atretic follicle (right). The horizontal axis 
represents the cell age, the vertical axis the cell maturity, and the colorbar indicates the cell density value 
(number of cells per elementary age x maturity volume). 02 is the cell cycle duration (2 age units). 

boundary value problems have been proposed (see e.g. [1] and [5] where renormalized entropy solutions are 
defined for this type of problem). An additional problem faced here consists in defining the trace of the 
solution on the outward boundary dQ^[ k , in order to solve the next problems Pi^+i- It is not clear whether 
this trace can be defined here for such weak solutions. As Vf £ (L co {Qi : k)) 2 and div(u/) £ L co {Qi^); an 
alternative is to use the results of [6], to define a unique solution <fi^ £ L 2 {Qi t k) with well defined trace in 
L 2 (dQ^ k ) as soon as the boundary data are in L 2 (dQ^ k ). A drawback of this alternative is that, in our 
case, it is not known whether those solutions coincide with classical weak solutions. 

2.3 Characteristic equations 

The model deals with controlled partial differential equations. Up to now, no convincing control strategy 
for such velocity controlled conservation laws with integro-differential control terms is available from the 
literature. To tackle the control problems, we focus on the actions of the control terms on the characteristics 
of the conservation laws. As the solution of the conservation laws is L 2 , we can derive the characteristic 
equations. 

We now deal with the following ordinary differential equations [7j: 

!dc =9f{lc,Uf) 
ic = h f ( lc ,u f ) c=l...n (20) 

0/ c = -(A( 7c ,C/) + ^4r £l ) 

where g/,hf and A arc defined in ([2]) , © , HJ) , (O and n is the number of characteristics. 

The local control Uf acts on the position of the characteristics on the (a c , j c ) spatial domain at a given time, 

and U on the loss term A. 

The transfer conditions between the cellular phases arc continuous on a c and 7 C and discontinuous on <pf c . 
From phase 1 to phase 2 (crossing the right boundary of some fii.fc), the flux continuity condition reads: 
- for T k such as: a c (T k ) = ax + ka 2 , < 7 c (? 1 i ) < 7 S , for some k £ N, 

0/cC#) = (l-3i^_(7)(l-M^7)))MT-) (21) 

From phase 2 to phase 1 (crossing the left boundary of some Oi,fe), after mitosis, the flux doubling condition 
reads: 
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- for Tk such as: a c (Tk) = ka^, < j c (Ti) < J s for some k £ N, 

(l- fllW ._(7)(l-u / (T+)))0 /c (2]+) = 20 /c (T-) (22) 

Ovulatory follicles are distinguished from the atretic ones by their characteristic curves reaching a well 
defined zone of the (a c ,7 c ,0/ c ) space at a given time, as we can see on Figure 2] In the next section, we 
study to which extent characteristics curves can reach the target sets when subject to appropriate open-loop 
control laws. 

3 Backwards reachable sets 

3.1 Link with Hamilton- Jacobi-Bellman equations 

We assume that a follicle ovulates (respectively undergoes atresia) at time T if: 

Vc = 1 . . . n, (a c (T), 7c(T), 4>f c (T)) £ M ( resp. £ M a ) where M and M a are the targets for respectively 
ovulation and atresia. 

We focus here on the backwards problem, and we first introduce the elements of its solution, the solvability 
tube [5] (or the backward reachable set [DJ) associated with M. D (resp. M a ), i.e. the states from which it is 
possible to reach M. Q (resp. Ai a ), given admissible controls U and uj. From a physiological viewpoint, it 
corresponds to the initial conditions compatible with ovulation (resp. atresia) and subject to correct control 
laws amongst the set of admissible controls. 

We recall briefly how such backwards reachability problems can be solved in the framework of optimal 
control theory. Let us consider the following nonlinear continuous controlled dynamics, where / is a Lipschitz 
continuous function, so that there is a unique continuous solution for any measurable u(t): 

X = f(t, X,u), t > T 

x{t) £ E" (23) 
u £ U C M" 1 

Given a closed target set M E R n and two times r and ti, the solvability set W(t, t\,M) is the set of states 
x G R n such that there exist control functions u with values in the compact set U that steer system (|23p from 
the state x(t) = x to x(t\) S M. [8]. W(t, t\,M) can be computed by solving the following optimization 
problem where d 2 (x,X) = min{||a; — z\\ 2 \ z £ X} is the square of the distance d(x,X) from point x to set 
X: 

min{d 2 (x(£i),.M)|2;(T) = x} (24) 

The solution of this optimal control problem can be found by solving the following backward Hamilton- 
Jacobi-Bellman (HJB) equations, where Hf(D x V,x) = min{(D x V) T .f(t, x,u)} is the Hamiltonian [ID] : 

V t +H f (D x V,x) =0, t<t u V{t l ,x)=d 2 (x{t l ),M) (25) 

Then the following property is true [8] : 

W(r,ti,M) = {x\V(t,x) < 0} (26) 

The solvability tube is then the set- valued function of the starting time t for given t\ and A4: 

W[h,M}: t^W[t 1 ,M](t) = W(t,t 1 ,M), T<t<h. (27) 

Any initial time t is thus associated with a delay of exactly (ti — t) to reach the target. 

In contrast, the backward reachable set is the set of states from which it is possible to reach the target in at 

most (ti - t) [DJ: 

G(t;ti,M)= |J W(t,s,M), t<h. (28) 

t<s<ti 
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This set can also be characterized in the same manner as in expression (f2"S|) . It is the zero sublevel set of 
the V function solution of the HJB equation: 



V t + mm(0,Hf(D x V,x)) 
Q(t ;h t M) 



0, t<h, V(h,x) = d 2 {x{ti),M) 
{x e R n \V{t,x) < 0} 



(29) 
(30) 



In the next section, we use the numerical methods based on equations (|2"5|) . (|3H)) to compute the Q set [TTlll2j . 
This choice is convenient since it seems easier to compute the Q set than the W set- valued function, but it 
amounts to deal with desynchronized final times. Indeed, in the W[ii,A4] tube, every states considered at 
the same starting time reach the M. target at the same final time t\, whereas in Q(t ;t±, M.) states reaching 
the M. target at different final times (provided that they arc not later than t\) arc fused. To ensure that 
every states reaching the target remain inside until t\, the control function has to be extended. 

3.2 Application to granulosa cells 

3.2.1 The follicle as a controlled dynamical system 

We represent each follicle by n cells following characteristic curves of (I} defined by three state variables: 
age, maturity and cell density as in ([20]) . The model of the follicle is then of the form (|23|) with 3ra state 
variables: 

x(t) = (oi(t), ...,a n (i);7i(*), ...,7„(i);0/i(i), ...,cf)f n {t)) T 
and for each cellular phase i = 1,2,3, the dynamics is defined by the following ff. 



( 5/(7i ,«/) \ 

9f{ln,Uf) 

h f (11,11 f) 

hf(i„,u f ) 
K(n,u f ,U)(j)f 1 



\ K(y n ,u f ,U)<f> fn j 
where <?/ and hf are defined in @, ([3]), ([5]) and K(~f c ,Uf,U) is defined by 



Phases 1 and 3: 



Phase 2: 



(31) 



Vfc, TV £ N, V(oc, 7c) G Qi, k U n 3 ,N 
K(7 C) , /)f /) = -(A(7 C) C/) + ^^ 
Vfc £ N, V(a c ,7 c ) G Q, 2 ,k 

Mk+l)a 2 

K(-f c ,Uf,U) =\n(2)5(a c )(f)f c , with / S(a)da = 1 



(32) 



(33) 



The expression of K(j c , Uf, U) in phase 2 results from the regularization by a smooth positive function 5(a) 
of the mitosis at the transition (|22p between phases 2 and 1, that will be used in the numerical treatment 
of the equations where the transmission boundaries are included in the computational domain. The precise 
choice for 6(a) does not matter since its only role is to guarantee that the density of cells entering phase 2 
at time ka\ has doubled at time ka\ + Since a 2 — <Xi = 1, the simplest choice is 5(a) = l so that we have 
6(a c ) = l, resulting in K(j c ,Uf,U) = ln(2)</> /c . 

According to physiological knowledge, wc call consider thctt there exists ci ma-ximal cige, o>maxi 

beyond which 

all cells become senescent and die. This assumption introduces the notion of a maximal life time T max such 
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that a(T max ) = a rnax (T max exists because a > g 2 > 0). Hence, there is also a maximal maturity 7 maa; and 
a maximal cell density 4>f m ax attainable for t G [0,T TOa J. The are Lipschitz continuous, so that there is 
a unique solution x within each phase. 

We now state some properties of this controlled system and the corresponding reachable targets in the 
age-maturity plane. We first describe the control action on the orientation of the velocity field: 

In the interior of phases 1 & 3: Vu/ E U, 

9h > giT g f = 0.5r g/ (34) 
hf("f,Uf) > 0, for Uf > u*j and < 7 < 7+(l/u) 

hf("f,Uf) < 0, for Uf < u*j and < 7 < 7 + (l/u) (35) 
hf{l,u f )<0, for all u f and j + (U/u) < 7 < ~f ma x 

In the interior of phase 2: Vit/ € U, 

g h = l, h f = 0. (36) 

with: 

*/ \ u\n ( 1 ^ 2 — rr ) for < 7 < 7+(l/?2) 

[ 1 for 7+(l/w) < 7 < 7mai 

7 ± M = C1{V) ± VCli f + ^ , with „(„) = c.(l - ex P (-,)), i=l,2. (38) 

We just sketch the proof (easy but tedious) of these properties. We can first remark that, in the interior of 
phases 1 and 3 excluding the strip defined by 7 S± , the u s± functions can only take as values or 1. In this 
interior subset hf{^f,Uf) is defined as 

hf(~f,Uf) = T hf (-"f 2 + (C17 + c 2 )(l - w s (7)exp(-tt//u))) 
= T h f(j+(uf/u) -j)(j-7-{uf/u)) 

It is easy to check that hf has the sign of (7 + (u//it) — 7). Besides, 7+ is an increasing function of it/, since 

, / x C ll+ {v) + C 2 

?+(") = y p^l = a =Ft exp(-i/) 
VciW^ + 4c 2 (^) 

Hence < 7+(u//u) < 7+(u m aa;/u) = 7+(l/?2), so that we can solve /i/(7*,u/) = and define 1*^(7) such 
that 7 + (mJ(7)/u) = 7 for < 7 < 7 + (1/u). Finally, as exp ( — — J is very small, 



lt/-\ ft I t l~\\ Cl + V C 1 + 4c 2 ^ 

7 + (1/m) « (1 - cxp(-l/u)) ^- < -y max 

This yields property (f3"5)h The remaining properties (|3"4"|) and (f3"6"| are obvious. 
Such results imply the following properties of the reachable targets: 



1. As 5^ > 0, the ages of cells in the reachable targets set have to be larger than those of the cells in 
their initial conditions. Hence, the reachable targets are on the right of the cell initial positions in the 
age-maturity plane. 

2. For a given global control U and maturities < 7 < j+(U/u), property (|33|) shows that it is possible, 
using a local feedback Uf around the feedback u*p to either increase, decrease or maintain the maturity. 
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In particular, if a target maturity is reached and U remains large enough, the cells remains inside the 
target. For such targets M, Qit ;t\,M.) — W(t ;t\,M.). It is also possible to reach a target and then 
leave it, e.g. if U is lowered. This set-equality depends then strongly on the control strategies and 
positions of the targets with respect to 7 + (l/u). 

3. As U < 1, maturities y+(l/u) < 7 < ^ m ax are not reachable. 
Subsequently, we assume that: 



so that it is not impossible to reach them (but success is not guaranteed, the density conditions needing also 
to be met). 

3.2.2 The Hamiltonian in each cellular phase 

The HamiltonJacobi-Bcllman (HJB) equations associated with the backwards reachable sets are solved once 
a target has been defined and the Hamiltonians in each phase have been minimized according to the control 
terms. We do not use here the feedforward term: Uo = 0. Hence both the global and local FSH signals 
operate as best-case feedback controls subject to the explicit constraints (rewriting (fT5)0 < Uf < U < 1. 
We note u = (uf, U) T the control and U the set of admissible controls. Such constraints take into account 
the following physiological knowledge: 

(i) the level of FSH in the antral fluid is at much as high as FSH plasmatic levels, which implies Uf <U: 

(ii) there is a basal secretion (i.e. not subject to ovarian feedback) of FSH by the pituitary gland, which 
implies U > 0; 

(iii) the balance between the secretion rate (or exogenous administration rate) and the clearance rate ensure 
the saturation of FSH plasmatic levels, which implies U < 1. 

Two coupled biological questions underlie this reachability problem: (1) which patterns of exogenous FSH 
administration can be applied to target either ovulation or atresia and (2) how follicular vascularization 
should develop to remain compatible with those patterns. The answer to such questions might help improving 
the current FSH treatments, which suffer from several drawbacks such as ovarian overstimulation syndrome. 

To reach a given target M. from its current location in phase 0,^ 0, a cell (or the lineage of its daughter 
cells) has to reach some intermediate or final target in the closure of the fij^ domain. 

For instance, in the case of the ovulatory target M = M {M C ^3,jv) W[*i, M ] following the solvability 
tube backward from Phase 3 leads to the boundary between and fli.k (k < N); the intersection of this 
boundary with the solvability tube defines a new target set M. ,i,k, in ^i,fc- 

Such intersections depend smoothly upon the width of the thin strip < "f s _ < 7 < 7 S+ containing these 
boundaries where the control has little effects. In the following, we suppose that 7 S _ = 7 = j s+ so that uj s _ 
and lo s+ arc replaced by perfect switch functions. 

In each cellular phase i, in order to compute the backward reachable set Q(ti — t ;ii,_M;), i.e. the set of 
states from which it is possible to reach the target Aii in less than a given time t, we have to solve an HJB 
equation of the form ([29]) and use the definition (f30|) : 



the target ages are large enough (e.g. > 01 + a 2 ) 
the target maturities are below 7 + (1/m) 



(39) 
(40) 



V t + min(0,H f i(D x V,x)) = 0, t < t u V(h,x) = d 2 (x(h),Mi) 
G(t;t!,M) = {x G R n \V(t, x) < 0} 



(41) 
(42) 



where Hfi(p, x) 



mm(p T .f t (x,u)), with p = D x V = (p ai , ...,p Qji ,p 7l , ...,p 7n ,p 0/1 , ...,p 0/ J T . 



u 



1 ovarian follicles are spheroidal structures hollowed by a cavity called the antrum 

2 where i stand for the phase index and k for the number of cell generation elapsed since initial time 
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If the targets fulfil the conditions and (T4T))) discussed above, we can conclude that: 

i) Starting from cell ages "younger" than the target ages, for all u €U, the cells get closer in age to target: 
For all i, p CH .9/(74, Uf) < 0, as p ai < and <?/ > 0. 

ii) For the particular choice U = 1 and Uf = u%, we have hf(-fi,Uf) = 0, so that: 
For all i, p 7i ./i/(7i, u*^) = 0. 

In particular, p T .fi(x, u*^ )) < 0, so that Hfi(p, x) < 0. We have thus shown that in our case, solving equation 
(|4"Tj) amounts to solve the following simpler HJB equation to compute Q: 

V t +H fi (D x V,x) = 0, t<h, V(t 1 ,x) = d 2 (x(t 1 ),M i ) (43) 
g{t;t x ,M) = {x G R n \V(t,x) < 0} (44) 

We now compute the Hamiltonian Hfi(p,x) in each phase. 
Phases 1 & 3 

We first compute p T -fi(x, u), for i = 1,3. 

n n n 

p T .f l {x 1 u) = ^Pa c 9f{u s ) + J2p lc hf{j c , Uf) + ^p0 /e 0/ c if (7 C , u f , U) 

c—1 c—1 c— 1 

which can be rewritten as follows: 

p T .f i (x,u) = H (p,x) - A(p,x)exp(-Uf/u) + B(p)u f + C(p,x)U 
with (using equations ([2]), ^ and ([52")) ') 

n n n 

H (p,x) = T gf (l - gi)^2pa c + T hf^2p lc {-Jc + C llc + C 2 ) + ^2p <j)fc (t)f c K (j c ), 

c—1 c—1 c—1 

where #0(7) = -(^(7) + T hf (-2j + a)), 

n 

A(p, x) = T hf (Pia( c Uc + ca) - ciP4, fc <f)f c ) , 

c=l 

n n 

B (P) = T gf9l^2Pa c , C(p,x) = ^p 0/c /c W A (7 c ). 
c=l c=l 

Now we minimize this expression with respect to u G IA. 

Hti(p,x) = min(p T .fi(x,u)) = H (p,x) + min (B(p)uf - A(p, x) exp(~Uf/u) +C(p,x)U). 

u£U 0<Uf<U<l 

It is equivalent to minimize first over U with Uf < U < 1 and then over Uf with < Uf < 1 (we drop the 
arguments x,p for simplicity): 

Hfi~H n + min min (Bit*- — Acxp(— Uf/u) + CU) I 

0<Uf<l \Uf<U<l J 

When C > 0, we have H fi = H + min {{B + C)u f - Aexp(-Uf /u)). 

0<uf<l 

When C < 0, we have Hfi = Ha + C + min (Buf — Acxp(— Uf/uj). 

Q<uf<l 

We have then to solve min [Buf — Aexp(—Uf /u) ) with B = B + C| + . 

When A > (rcsp. A < 0), B«/ — Acxp(— Uf/u) is a concave (resp. convex) function of u/ and the 
minimizing control is now easy to compute, leading to the following feedback law: 

When A(p,x) > 0: u f = if - A(p,x) < B(p,x) - A(p, x) exp(-l/u), else u f = 1 (45) 
When A(p,x) < 0: u f = -u\n(- B ^ ,X '™ ) if f ^' ^ > exp(-l/u), else U/ = 1 (46) 
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Phase 2 

This phase is not controlled so that the expression of the Hamiltonian is straightforward using (|33[) : 

n n 

H f2 = T gf ^p ac +\YL{2)^2 P4>j .J{a c )(t>f c 

c=l c=l 

This analytical minimization enables to write the expression of the Hamiltonian in each cellular phase. 

3.3 Numerical results 

3.3.1 Implementation: level set methods 

Two main approaches are used to find reachable sets: the Eulcrian approach approximates the solution 
values on a fixed grid, whereas the Lagrangian approach tracks the trajectories of the dynamics. In case of 
backwards reachable sets, Eulerian methods are the most appropriate, as they enable to correctly compute 
the solution beyond shocks [13] . 

Amongst the various Eulerian methods available to compute the exact reachable sets, we used "level set 
methods" to solve Eq. (|25|) . They enable to simulate the motion of dynamic surfaces, such as the zero-level 
set of HJB equations |12j , with a very high accuracy (about a tenth of the spacing between the grid points 

m n 

Such level set algorithms are implemented in the "Toolbox of Level Set Methods" [j, which can be used to 
solve, among others, equations of the form: 

D t V(x,t) + H(x,D x V)=0 
with the constraints DtV(x, t) > or V(x, t) > V(x, t) 

The time derivative, D t V, is approximated by a Runge-Kutta scheme, with customizable accuracy. The 
spatial derivative, D X V, is approximated with an upwind finite difference scheme. The Hamiltonian, H(x,p), 
is approximated with a Lax-Friedrichs scheme: 

H(x,p + ,p~) = H(x, P ~l P ) - ^a T (p + -p~) 

where p + and p~ are respectively the right and left approximations of p, and a is an artificial dissipation 
coefficient added to avoid oscillations in the solution. 

Although the toolbox is designed for solving initial value problems, it is possible, in case of autonomous 
systems, to solve final value problems by reversing the time in system (f23| p~2] . 

3.3.2 Simulation results 

The initial problem is a 3 x n dimensional problem, where n is the number of characteristics considered. 
This problem is not tractable from a numerical ground. Indeed, for n > 2 and a correct number of grid 
points in each dimension (~ 100) we cope with too many points (100 3 "), and face the dimensionality limits 
of the algorithms. We thus computed the backwards reachable sets for one characteristic, (a±, 71, </>/i), 
corresponding to an ovulatory or atretic trajectory. It is a simplified case where all granulosa cells in a 
follicle are assumed to be synchronized in age and maturity. 

In the case of ovulatory follicles, the boundaries of the target set M are chosen to coincide with the ranges 
of age, maturity, and cell density reached by the characteristics of an ovulatory follicle at ovulation time, 
(see Figure |H left). 

The target set A4 for ovulation is thus defined as: 

Mo = {(01,71,^/1) S [10,12] x [10,11] x [4,6]} (47) 

3 http: / /www.cs.ubc.ca/~mitchell/ToolboxLS/ 
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Figure [5] illustrates the trajectories compatible with ovulation. 

The left panel illustrates the projection, on the age-maturity plane, of the changes in the backwards reachable 
set. Time t = (first subplot) is ovulation time, and we can see the rectangle {[10, 12] x [10, 11]} representing 
the projection of the target set M. (the age range is slightly different from that of Figurc[l]duc to difference 
in the dimensioning of the unrolled domain compared to the periodic domain). At time t = 4, the set area 
has increased, including all states that can join the target set in less than 4 time units (2 time units roughly 
correspond to 1 cell cycle duration). The backwards reachable set at time t = 11 represents the set of states 
that can reach the target set for ovulation in at most 11 time units. We chose to stop the simulation at this 
time as it contains the initial conditions used for the simulation represented in Figure |H The "staircase" 
shape of the backwards reachable set is due to the dynamics in the cell cycle. The maturity increases with 
age in phase 1, while its remains unchanged in phase 2, yielding to the flat part of the staircase. An instance 
of the flat part can be seen at time t = 11 between a\ = 7 and a\ = 8. 

The right panel illustrates the projection, on the age-cell density plane, of the changes in the backwards 
reachable set. We can again see the projection of the target set M at time t = 0, as the rectangle 
{[10,11] x [4,6]}. At time t = 11, a large part of the age-cell density domain is compatible with reaching 
the target set. From age <xi = 9 onward, the cell density variable increases with age to reach the target set. 
This is in agreement with the cell density dynamics in phase 3. 

The backwards reachable set of M Q nearly contains the whole spatial domain. It is an overapproximation of 
the initial states that allow to reach the target set. Indeed, it does not only contain the trajectories of the 
characteristics of the ovulatory follicle, but also those starting from initial conditions that do not arise on a 
physiological ground. Basing on biological knowledge, we can thus select from the overapproximated set, the 
physiologically-meaningful subset: we assume that all cells arc initially within the cell cycle, at their first 
generation, so that the subset of admissible initial conditions is defined by {(ai(io), 7i(^o)) £ [0, a2\ x [0, 7s]}- 
This subset is below the dashed segment on the left panel of Figure [5] 




Figure 5: Solvability tubes for ovulatory follicles. The left panel shows the projection of the 3D set on the 
(ai-71) plane, and the right one on the (ai-0/i) plane. 

The target set M. a for atresia is defined by the ranges reached by the characteristics of atretic follicles: 

M a = {(01,71,^/1) G [8,10] x [3,4] x [2,4]} (48) 

Figure [5] illustrates all the trajectories that fall into the target set for an atretic follicle. The left panel 
illustrates the projection, on the age-maturity plane, of the changes in the backwards reachable set. Some 
states in the target set are initially in the cell cycle, leading to the "staircase" pattern described below. The 
right panel illustrates the projection, on the age-cell density plane, of the temporal changes in the backwards 
reachable set. The set area increases, and at time t = 11, it roughly contains the whole domain. Once again, 
the set of initial states leading to atresia is overapproximated, and we can apply the same constraints as in 
the ovulatory case on the initial conditions, that are again below the dashed segment on Figure [51 
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Figure 6: Solvability tubes for atretic follicles. The left panel shows the projection of the 3D set on the 
(01-71) plane, and the right one on the (ai-(/>/i) plane. 



The initial conditions leading to ovulation or atresia, subject to physiological constraints, are roughly su- 
perimposed, and contain the whole subset of admissible initial conditions. Thus we cannot distinguish a 
priori ovulatory trajectories from atretic ones. This result agrees with physiological knowledge, as there is 
no predestination in follicular fate [H] , 

Although the initial conditions do not enable to determine whether a follicle is ovulatory or atretic, there 
is a moment when trajectories leading to either ovulation or atresia diverge. Such a moment corresponds 
to the physiological time of selection, when the trajectories of the few ovulatory follicles separate from the 
atretic ones. 

4 Discussion and perspectives 

We have characterized the backwards reachable sets for both ovulatory and atretic follicles. This control 
problem has been solved by studying the characteristics of the conservation laws describing the dynamics of 
the structured granulosa cell populations in ovarian follicles. The HJB equations representing the backwards 
reachable sets have been solved numerically for one characteristic curve, and allow to characterize the set of 
initial conditions that lead a characteristic into the target set either for ovulation or atresia. 

We have used the backwards reachable sets theory for a problem that is not entirely continuous, as there 
are discontinuities at the transitions between each cellular phase. Hence, we studied each phase separately, 
assuming that the transitions at the boundaries were continuous. Actually, we deal with a hybrid problem, 
as each cellular phase has its proper dynamics. As far as time dependent HJB equations are concerned, 
elements of a theory for hybrid systems are established, especially concerning "reach-avoid" sets [15], but it 
is not complete yet, and not directly applicable to our problem. Viability theory also considers reachable 
sets for hybrid systems [16], but the algorithms used are not as accurate as level set methods, since they use 
a different representation of reachable sets [T3] . Yet the numerical results we have obtained with the level 
set methods seem correct since the continuous part of the system is solved with highly accurate algorithms 
and the hybrid part is handled via continuity conditions and substantiated by the trace results obtained on 
the conservation laws. 

We have solved a simplified numerical problem dealing with a single characteristic of a follicle, as the 
simulation tool needs too huge memory to solve higher dimensional problems. The backwards reachable sets 
obtained show that it is possible to find a correct control law to steer a large set of initial conditions into 
the target set, and they give information upon the maximum duration needed to reach the target. 
If the simulation of more than one characteristic were possible, we expect that the backwards reachable sets 
would be smaller. We even speculate that, from a physiological viewpoint, a larger set in case of a single 
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growing follicle ensures the ovulation success, while a smaller set in case of several simultaneously growing 
follicles limits the ovulation rate and underlies the selection process in part. Yet this assumption needs 
improvement of the numerical implementation to be tested. 

We have studied by now the open-loop problem, for the ovulation of one follicle. This situation is close to 
therapeutic situations of ovarian stimulations, so that the designed control laws may be useful in improving 
therapeutic schemes. Yet, to understand ovulation control in depth, we will have to tackle multi-scale, 
closed-loop reachability problems of selection, which is the matter of future work. 
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